// GPU Verlet Integrator
// Autor: Juan G. Campa

#ifndef _VERLETINTEGRATOR_CU_
#define _VERLETINTEGRATOR_CU_

// includes, C string library
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

// includes, cuda
#include <cuda.h>
#include <builtin_types.h>
#include <cuda_runtime_api.h>

#include <Common.cu>
// includes, project
//#include <cutil_inline.h>

#define clampValue(x,y,z) x<y?y:x>z?z:x
#define SPHERE_RAD 0.7

extern "C" float4 __device__  factor(float4 pos, float factor);

extern "C" __global__ void Integrate(float4 *posArray, float4 *prevPosArray, const int particleCount, const float dt)
{
   unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;

   /*float4 aux = posArray[i];
   float4 a = aux + aux;
   prevPosArray[i] = aux;*/
   
   if (i < particleCount)
   {
      float4 curPos = posArray[i];
      float4 temp = curPos;
      float4 prevPos = prevPosArray[i];
      float4 forces;
      forces.y = -0.0002f;
      forces.x = 0;
      forces.z = 0;
      forces.w = 0;

      //curPos = 2* curPos - prevPos + forces ;
      curPos = factor(curPos, 1.99) - factor(prevPos, 0.99) + forces;

      curPos.y = clampValue(curPos.y, -2, 2);
      curPos.x = clampValue(curPos.x, -3, 3);

      float distToCenter = sqrt(curPos.x * curPos.x + curPos.y * curPos.y + curPos.z * curPos.z);
      if (distToCenter < SPHERE_RAD)
      {
          // Normalizamos y multiplicamos por el radio de la esfera.
          curPos.x = curPos.x / distToCenter * SPHERE_RAD;
          curPos.y = curPos.y / distToCenter * SPHERE_RAD;
          curPos.z = curPos.z / distToCenter * SPHERE_RAD;
      }

      posArray[i] = curPos;
      prevPosArray[i] = temp;
   }
}

/// Multiplies the float4 by the given factor.
extern "C" 
__device__ float4 factor(float4 pos, float factor)
{
    pos.x = pos.x * factor;
    pos.y = pos.y * factor;
    pos.z = pos.z * factor;
    pos.w = pos.w * factor;
    return pos;
}

extern "C" 
void ExecuteIntegrate(void* pos, void* prevPos, unsigned int particleCount, float dt)
{
    cudaError_t error = cudaSuccess;

    dim3 block(256, 1, 1);
    int gridX = particleCount / block.x;
    if ((particleCount % block.x) > 0)
        gridX++;
    dim3 grid(gridX, 1, 1);

	Integrate<<<grid, block>>>((float4 *)pos, (float4 *)prevPos, particleCount, dt);

    error = cudaGetLastError();
    if (error != cudaSuccess) {
        printf("Integrate() failed to launch error = %d\n", error);
    }
}
#endif